Baseline Cumulative DS

• Cumulative David’s Score based on chases

res = lapply(1:4, function(d) {
  p = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:d), wins = 'Chases', nPerm = 1,
                        idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
    dplyr::select(GroupNumber, Sex, DS) %>% unnest() %>% 
    mutate(MouseID = as.numeric(MouseID)) %>% rename(DS_cum = DS)
  p$Day = d
  return(p)
})
res = do.call('rbind', res)
res = left_join(res, obj$mice[,c('MouseID', 'GroupNumber', 'Sex', 'RealMouseNumber')])

# Rank day 4
p2 = res %>% filter(Day == 4) %>% group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-DS_cum), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )
res = left_join(res, p2[,c('RealMouseNumber', 'Rank')])

r2 = res %>% group_by(Sex, Day, Rank) %>% 
  dplyr::summarize(
    M = mean(DS_cum, na.rm = T),
    se = sd(DS_cum, na.rm = T)/sqrt(n()),
  ) 
f1_1 = ggplot(r2, aes(x = as.factor(Day), y = M, fill = Rank, 
                      color = Rank)) + 
  geom_line(data = res, aes(group = RealMouseNumber, y = DS_cum), alpha = 0.4, size = 0.4) + 
  geom_errorbar(aes(ymin = M-se, ymax = M+se), width = 0) +
  geom_line(aes(group = Rank), size = 0.8) + 
  labs(y = "Baseline David's Score", x = 'Day') +
  geom_point(size = 2.5, pch = 21, color = 'grey20') + 
  scale_fill_manual(values = obj$colors_rank, name = 'Final rank', 
                    labels = c('\u03B1', '\u03B2', '\u03B3', '\u03B4')) + 
  scale_color_manual(values = obj$colors_rank, name = 'Final rank', 
                     labels = c('\u03B1', '\u03B2', '\u03B3', '\u03B4')) + 
  scale_x_discrete(expand = c(0.1, 0.1)) +
  facet_grid(~Sex)  + theme(legend.position = 'top')
f1_1

Correlations with SB behaviors

• Cumulative David’s Scores from days 1:4
• SB behaviors summarized as means of days 2:4

wins = 'Chases'
h = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = wins, nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = F)
h = h %>% dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(as.character(MouseID)))
h = left_join(h, obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])

p = subset(obj$profiles, !Suspect & Day %in% 2:4) %>% 
  dplyr::select(RealMouseNumber, all_props) %>% 
  group_by(RealMouseNumber) %>% 
  dplyr::summarise_at(.vars = all_props, mean, na.rm = T)

d = left_join(h, p)

idVars = c('RealMouseNumber', 'GroupNumber', 'Sex', 'DS')
props = c('DistanceOutside', 'GridEntropy6x6', 'FractionOfTimeOutside')

p <- d  %>% 
  dplyr::select(idVars, props) %>% 
  gather(prop, value, -idVars) %>% 
  mutate(prop = plyr::mapvalues(prop, all_props, all_props_pretty))

f1_2 = ggplot(p, aes(x = DS, y = value, fill = Sex)) +
  geom_smooth(method = 'lm', se = F, lty = 2, aes(col = Sex)) +
  geom_point(pch = 21, size = 2) +
  labs(x = "Baseline David's Score", y = '') +
  facet_grid(prop~Sex, scales = 'free_y') +
  scale_fill_manual(values = obj$colors_sex) +
  scale_color_manual(values = obj$colors_sex) +
  theme(legend.position = 'none')
f1_2

wins = 'Chases'
h = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = wins, nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = F)
h = h %>% dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(as.character(MouseID)))
h = left_join(h, obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])

p = subset(obj$profiles, !Suspect & Day %in% 2:4) %>% 
  dplyr::select(RealMouseNumber, all_props) %>% 
  group_by(RealMouseNumber) %>% 
  dplyr::summarise_at(.vars = all_props, mean, na.rm = T)

d = left_join(h, p)

idVars = c('RealMouseNumber', 'GroupNumber', 'Sex', 'DS')
props = all_props

p <- d  %>% 
  dplyr::select(idVars, props) %>% 
  gather(prop, value, -idVars) %>% 
  group_by(prop, Sex) %>% mutate(N = n()) %>%
  group_by(prop, Sex, N) %>% nest() %>% 
  mutate(
    M = map(data, function(d) {
      tidy(cor.test(d$value, d$DS, method = 'spearman'))
    })
  ) %>% dplyr::select(-data) %>% unnest() %>% 
  mutate(
    prop = fct_reorder(prop, abs(estimate))
  ) %>% group_by(Sex) %>% 
  mutate(qval = p.adjust(p.value, method = 'BH')) %>% 
  arrange(p.value) %>% 
  mutate(sig_qval = qval <= 0.1) 

knitr::kable(table(p$prop[p$sig_qval], p$Sex[p$sig_qval]))

p[p$prop == 'DistanceOutside',]

p[p$prop == 'FractionOfTimeOutside',]

p[p$prop == 'GridEntropy6x6',]

DS time correlations

• Pearson correlations of DS values between Day 1 and the following days

p = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = T) %>% 
  dplyr::select(GroupNumber, Sex, Condition, Day, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(MouseID)) %>% 
  group_by(GroupNumber, Sex, Condition, Day) %>% 
  mutate(DS_rank = rank(-DS))

p = left_join(p, obj$mice[,c('RealMouseNumber', 'GroupNumber', 'Sex', 'MouseID')])
p$Stage = ifelse(p$Day %in% 1:4, 'Baseline', 'Stress')

p = p %>% group_by(GroupNumber, Sex, Condition, MouseID) %>% arrange(GroupNumber, Sex, MouseID, Day) %>% 
  dplyr::select(-DS) %>% rename(DS = DS_rank) %>% 
  mutate(
    lag0 = DS, lag1 = lag(DS), lag2 = lag(DS, 2), lag3 = lag(DS, 3), lag4 = lag(DS, 4), 
    lead0 = DS, lead1 = lead(DS), lead2 = lead(DS, 2), lead3 = lead(DS, 3), lead4 = lead(DS, 4)
  )
pp = p %>% group_by(Day, GroupNumber, Sex) %>% nest() %>% 
  mutate(
    lead = map(data, function(d) {
      x = cor(d[, c('lead0', 'lead1', 'lead2', 'lead3')], 
              use = "pairwise.complete.obs", method = 'pearson') %>% 
        melt() %>% filter(Var1 == 'lead0')
    })
  ) %>% dplyr::select(-data) %>% unnest() %>% 
  mutate(uGroup = paste0(GroupNumber, '_', Sex))

pd = 0.2
pp %>% group_by(Var2, Sex, Day) %>% 
  dplyr::summarize(
    M = mean(value, na.rm = T), 
    se = sd(value, na.rm = T)/sqrt(n())) %>%
  filter(Day == 1) %>% ungroup() %>%
  mutate(Var2 = plyr::mapvalues(Var2, c('lead0', 'lead1', 'lead2', 'lead3'), c(1:4))) %>%
  ggplot(aes(x = Var2, y = M, color = Sex, fill = Sex)) + 
  # geom_hline(yintercept = 0, lty = 2, color = 'grey') +
  geom_line(aes(group = Sex), position = position_dodge(pd)) +
  ylab("David's Score\ncorrelation with day 1") + xlab('Day') + facet_grid(~Day, scales = 'free_x') +
  scale_fill_manual(values = obj$colors_sex) +
  scale_color_manual(values = obj$colors_sex) +
  geom_errorbar(aes(ymax = M+se, ymin = M-se), width = 0, position = position_dodge(pd)) +
  geom_point(size = 2, pch = 21, color = 'grey20', position = position_dodge(pd))  +
  theme(legend.position = 'top')

Sum of daily changes

# Generate a distribution of expected rank changes
set.seed(42)
sampleGroupRanks = function (i) {sample(1:i, size = i, replace = F)}

# Fake data
nMice = length(unique(obj$master$RealMouseNumber))
nDays = 4
nIter = 10

ref = lapply(1:nIter, function(i) {
  dat = data.frame(
    Day = rep(1:nDays, nMice),
    MouseNumber = rep(1:nMice, each = nDays), 
    GroupNumber = rep(1:(nMice/4), each = 4*nDays)
  )
  dat = dat %>% group_by(GroupNumber, Day) %>% 
    nest() %>%
    mutate(
      Rank = map(data, function(d) sampleGroupRanks(nrow(d)))
    ) %>% unnest()
  
  d = dat %>% group_by(GroupNumber, MouseNumber) %>% 
    arrange(GroupNumber, MouseNumber, Day) %>% 
    mutate(
      lagRank = lag(Rank, 1),
      RankChange = abs(Rank - lagRank)
    ) %>% filter(!is.na(lagRank)) 
  ref = d %>% group_by(GroupNumber, MouseNumber, RankChange) %>%
    tally() %>%
    group_by(GroupNumber, MouseNumber) %>% 
    dplyr::mutate(
      s = sum(n),
      SumRankChangesFraction = n / sum(n) 
    ) %>% 
    group_by(RankChange) %>% 
    tally() %>% 
    dplyr::mutate(
      SumRankChangesFraction = n / sum(n)
    )
  ref$iter = i
  return(ref)
})

ref = do.call('rbind', ref)
ref = ref %>% group_by(RankChange) %>% 
  dplyr::summarise(
    M = mean(SumRankChangesFraction, na.rm = T), 
    Qt = qt(0.975,n()-1) * sd(SumRankChangesFraction)/sqrt(n()),
    CI_95 = M + Qt,
    CI_05 = M - Qt
  )

# Real final ranks
r = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(MouseID)) %>% rename(DS_cum = DS)
r = left_join(r, obj$mice[,c('MouseID', 'GroupNumber', 'Sex', 'RealMouseNumber')])

r = r %>% group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-DS_cum), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )

# Real ranks throughout
p = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = T) %>% 
  dplyr::select(GroupNumber, Sex, Condition, Day, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(MouseID)) %>% 
  group_by(GroupNumber, Sex, Condition, Day) %>% 
  mutate(DS_rank = rank(-DS))
p = left_join(p, obj$mice[,c('RealMouseNumber', 'GroupNumber', 'Sex', 'MouseID')])
p$Stage = ifelse(p$Day %in% 1:4, 'Baseline', 'Stress')
p = left_join(p, r[,c('RealMouseNumber', 'Rank')])
p = p %>% group_by(GroupNumber, Sex, Condition, MouseID, Rank) %>% 
  arrange(GroupNumber, Sex, MouseID, Day) %>% 
  mutate(
    lagRank = lag(DS_rank),
    RankChange = abs(DS_rank - lagRank)
  ) %>% filter(!is.na(RankChange))

pp = p %>% 
  # group_by(Sex, Rank, RankChange, RealMouseNumber) %>% 
  group_by(Sex, Rank, RankChange) %>%
  tally() %>%
  # group_by(Sex, Rank, RealMouseNumber) %>% 
  group_by(Sex, Rank) %>%
  mutate(
    s = sum(n),
    SumRankChangesFraction = n / sum(n, na.rm = T) 
  ) 

ggplot(pp, aes(x = RankChange, y = SumRankChangesFraction)) +
  geom_ribbon(data = ref, aes(ymax = CI_95, ymin = CI_05, y = M), 
              fill = 'grey', alpha = 0.5) +
  geom_line(data = ref, color = 'grey', lty = 1, aes(y = M)) +
  geom_line(aes(col = Rank)) +
  geom_point(aes(fill = Rank), pch = 21, color = 'grey20', size = 2.5) +
  facet_grid(.~Sex) + 
  labs(y = 'Probability') +
  scale_color_manual(values = obj$colors_rank) +
  scale_fill_manual(values = obj$colors_rank) + 
  theme(legend.position = 'top')

x = left_join(pp, ref)
b = seq(-2, 2, 1)
l = c(-4, -2, 0, 2, 4)
ggplot(x, aes(x = as.factor(RankChange), y = log2(SumRankChangesFraction / M))) +
  geom_hline(yintercept = 0, lty = 2, color = 'grey') +
  geom_line(aes(col = Rank, group = Rank), size = 0.5, 
            position = position_dodge(0.4)) +
  geom_point(aes(fill = Rank, col = Rank), pch = 21,
             color = 'grey20',
             size = 2.5, 
             position = position_dodge(0.4)) +
  labs(y = 'Fold change in probability\n(true / expected)', x = 'Absolute rank change') +
  facet_grid(.~Sex) + 
  scale_y_continuous(breaks = b, labels = l) +
  scale_fill_manual(values = obj$colors_rank) + 
  scale_color_manual(values = obj$colors_rank) + 
  theme(legend.position = 'none')

Change/no change version

# Real final ranks
r = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(MouseID)) %>% rename(DS_cum = DS)
r = left_join(r, obj$mice[,c('MouseID', 'GroupNumber', 'Sex', 'RealMouseNumber')])

r = r %>% group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-DS_cum), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )

# Real ranks throughout
p = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = T) %>% 
  dplyr::select(GroupNumber, Sex, Condition, Day, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(MouseID)) %>% 
  group_by(GroupNumber, Sex, Condition, Day) %>% 
  mutate(DS_rank = rank(-DS))
p = left_join(p, obj$mice[,c('RealMouseNumber', 'GroupNumber', 'Sex', 'MouseID')])
p$Stage = ifelse(p$Day %in% 1:4, 'Baseline', 'Stress')
p = left_join(p, r[,c('RealMouseNumber', 'Rank')])

p = p %>% group_by(GroupNumber, Sex, Condition, MouseID, Rank) %>% 
  arrange(GroupNumber, Sex, MouseID, Day) %>% 
  mutate(
    lagRank = lag(DS_rank),
    RankChange = abs(DS_rank - lagRank)
  ) %>% filter(!is.na(RankChange))

pp = p %>% 
  dplyr::mutate(
    NoChange = RankChange == 0
  )

ppp = pp %>% group_by(Rank, Sex) %>%
  dplyr::summarise(
    N = n(), 
    S = sum(NoChange)
  ) %>%
  dplyr::mutate(
    P = S / N
  )

ppp$sig = NA
ppp$sig[ppp$Rank == 'alpha'] = '***'
ppp$sig[ppp$Rank == 'delta' & ppp$Sex == 'Male'] = '***'


ppp$Rank = plyr::mapvalues(ppp$Rank, c('alpha', 'beta', 'gamma', 'delta'), 
                           c( '\u03B1', '\u03B2', '\u03B3', '\u03B4'))
f1_3 = ggplot(ppp, aes(x = Rank, y = P / .25, fill = Rank)) +
  geom_bar(stat = 'identity', width = 0.7, alpha = 1, aes(color = Rank)) +
  geom_hline(yintercept = 1, lty = 1, color = 'grey30') +
  geom_text(aes(label = Rank, y = 1), nudge_y = 0.08, parse = F, fontface = 'italic', color = 'grey30') +
  geom_text(aes(label = N), color = 'grey20', nudge_y = 0.08, parse = F, size = 3, fontface = 'italic', color = 'grey30') +
  geom_text(aes(label = sig), color = 'grey20', nudge_y = 0.15, parse = F, size = 4, fontface = 'bold', color = 'grey30') +
  facet_grid(.~Sex) + 
  labs(y = 'Rank maintenance odds', x = '') +
  scale_color_manual(values = obj$colors_rank) +
  scale_fill_manual(values = obj$colors_rank) +
  scale_y_continuous(trans = 'log2', breaks = c(0.66, 1, 1.5, 2, 2.5, 3), limits = c(0.66, 3)) +
  theme(legend.position = 'none', 
        axis.text.x =  element_blank())

f1_3

Stats:

res = lapply(1:nrow(ppp), function(i) {
  d = ppp[i,]
  tidy(binom.test(d$S, n = d$N, p = .25, alternative = 'greater'))
})
res = do.call('rbind', res)
ppp = cbind(as.data.frame(ppp), as.data.frame(res)) %>% arrange(Sex, Rank)
as.data.frame(ppp)

DS baseline v stress

baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                             idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest()

stress = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 5), wins = 'Chases', nPerm = 1,
                           idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, DS) %>% unnest() %>% 
  rename(DS_stress = DS)

h = left_join(baseline, stress)

Baseline v. Stress DS

f1_4 = ggplot(h, aes(x = DS, y = DS_stress, fill = Sex)) +
  geom_smooth(method = 'lm', se = F, lty = 2, aes(col = Sex)) +
  geom_point(pch = 21, size = 2) +
  labs(y = "David's Score\nfollowing acute stress", x = "Baseline David's Score") +
  facet_grid(.~Sex) +
  scale_fill_manual(values = obj$colors_sex) +
  scale_color_manual(values = obj$colors_sex) +
  theme(legend.position = 'none')
f1_4

h %>% group_by(Sex) %>% 
  filter(!(is.na(DS) | is.na(DS_stress))) %>% 
  dplyr::summarize(n = n())

hh = subset(h, Sex == 'Male')
data.frame(tidy(cor.test(hh$DS, hh$DS_stress)), N = nrow(hh))

hh = subset(h, Sex == 'Female')
data.frame(tidy(cor.test(hh$DS, hh$DS_stress)), N = nrow(hh))

Chases after stress

baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 2:4), wins = 'Chases', nPerm = 1,
                             idVars = c('Sex', 'Condition'), verbose = F, daily = T) %>% 
  dplyr::select(GroupNumber, Day, Sex, Condition, data) %>% unnest() %>% 
  group_by(GroupNumber, Sex, Condition, MouseID, Day) %>% 
  dplyr::summarize(wins = sum(wins, na.rm = T)) %>% # Collapse over OtherID
  group_by(GroupNumber, Sex, Condition, MouseID) %>% 
  dplyr::summarize(wins = mean(wins, na.rm = T)) # Collapse over Day

stress = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 5), wins = 'Chases', nPerm = 1,
                           idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, data) %>% unnest() %>% 
  rename(wins_stress = wins) %>% 
  group_by(GroupNumber, Sex, MouseID) %>% 
  dplyr::summarize(wins_stress = sum(wins_stress, na.rm = T))  # Collapse over OtherID

h2 = left_join(baseline, stress)

Chases from baseline to stress

• Mean number of chases daily: baseline (days 2:4) vs. stress

pp = h2 %>% group_by(GroupNumber, MouseID, Sex, Condition) %>% 
  summarise_at(.vars = c('wins', 'wins_stress'), mean, na.rm = T) %>% 
  gather(Stress, Chases, -GroupNumber, -Sex, -Condition, -MouseID) %>% 
  mutate(Stress = plyr::mapvalues(Stress, c('wins', 'wins_stress'), c('Baseline', 'Acute Stress')), 
         Stress = factor(Stress, levels = c('Baseline', 'Acute Stress')))

Md = pp %>% group_by(Sex, Stress) %>% 
  dplyr::summarise(
    Md = median(Chases, na.rm = T),
    Mn = mean(Chases, na.rm = T),
    SE = sd(Chases, na.rm = T) / sqrt(n()),
    MAD = stats::mad(Chases, na.rm = T)
  )
# ppp = left_join(pp, Md)
f1_5 = ggplot(pp,
              aes(x = Chases, fill = Sex, color = Sex)) +
  geom_density(alpha = 0.7) +
  geom_rug(aes(color = Sex), size = 1, show.legend = F) +
  geom_errorbarh(data = Md[Md$Sex == 'Male',], inherit.aes = F, 
                 aes(xmin = Md, xmax = Md + MAD, y = -0.1, col = Sex), height = 0.1, show.legend = F) +
  geom_point(data = Md[Md$Sex == 'Male',], 
             pch = 21, size = 2.5, aes(x = Md, y = -0.1), color = 'grey20', show.legend = F) +
  geom_errorbarh(data = Md[Md$Sex == 'Female',], inherit.aes = F, 
                 aes(xmin = Md, xmax = Md + MAD, y = -0.2, col = Sex), height = 0.1, show.legend = F) +
  geom_point(data = Md[Md$Sex == 'Female',], 
             pch = 21, size = 2.5, aes(x = Md, y = -0.2), color = 'grey20', show.legend = F) +
  scale_fill_manual(values = obj$colors_sex) +
  scale_color_manual(values = obj$colors_sex) +
  facet_grid(Stress~.) + 
  scale_x_log10() + labs(y = '') +
  theme(legend.position = c(0.1, 0.9), 
        legend.title = element_blank(),
        panel.grid.major.x = element_line(color = 'grey', size = 0.2, linetype = 2),
        # panel.grid.minor.x = element_line(color = 'grey', size = 0.1),
        panel.border = element_blank(), 
        axis.text.y = element_blank())
f1_5

s = h2 %>% gather(Stage, Chases, -GroupNumber, -Sex, -Condition, -MouseID)
s$Stage = plyr::mapvalues(s$Stage, c('wins', 'wins_stress'), c('Baseline', 'Stress'))
s = left_join(s, obj$mice[,c('GroupNumber', 'Sex', 'MouseID', 'RealMouseNumber')])

summary(aov(log(Chases + 1) ~ Sex * Stage + Error(as.factor(RealMouseNumber)), data = s))

## 
## Error: as.factor(RealMouseNumber)
##           Df Sum Sq Mean Sq F value Pr(>F)
## Sex        1   0.34  0.3366   0.211  0.648
## Stage      1   0.59  0.5950   0.372  0.544
## Sex:Stage  1   0.01  0.0093   0.006  0.940
## Residuals 80 127.88  1.5984               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value      Pr(>F)    
## Stage      1  5.704   5.704  30.559 0.000000415 ***
## Sex:Stage  1  0.004   0.004   0.023       0.881    
## Residuals 78 14.559   0.187                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Figure panel

col1 = cowplot::plot_grid(plotlist = list(f1_1, f1_2), labels = c('a', 'b'), label_fontface = 'plain', 
                          ncol = 1, align = 'v', axis = 'lr', rel_heights = c(1, 2))
col2 = cowplot::plot_grid(plotlist = list(f1_3, f1_4, f1_5), labels = c('c', 'd', 'e'), label_fontface = 'plain', 
                          ncol = 1, align = 'v', axis = 'lr', rel_heights = c(1, 0.8, 1.2))
cowplot::plot_grid(plotlist = list(col1, col2), ncol = 2, align = 'h')

if(exportFigures) {
  ggsave('../manuscript/figures/1_Panel.pdf', width = 9.5, height = 9.3, dpi = 1200)
  ggsave('../manuscript/figures/1_Panel.png', width = 9.5, height = 9.3, dpi = 1200)
}

Physiological Outcomes

dat = obj$tests_combined_adj
dat$Sex = factor(dat$Sex, levels = c('Female', 'Male'))
dat$SexCondition = paste0(dat$Sex, '_', dat$Condition)
dat$SexCondition = factor(dat$SexCondition, 
                          levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))

dat$value = dat[['ORG_BodyweightChange']]
dat %>% group_by(Sex, Condition) %>% 
  dplyr::summarize(
    n = n(), 
    n_groups = length(unique(GroupNumber)),
    M = mean(value, na.rm = T), 
    sd = sd(value, na.rm = T)
  )

dat %>% group_by(Sex, Condition) %>% nest() %>% 
  mutate(SW = map(data, function(d) shapiro.test(d$value)), 
         tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()

knitr::kable(car::leveneTest(value ~ Sex * Condition, dat))

knitr::kable(car::leveneTest(value ~ Condition, dat[dat$Sex == 'Male',]))

knitr::kable(car::leveneTest(value ~ Condition, dat[dat$Sex == 'Female',]))

as.data.frame(tidy(aov(value ~ Condition * Sex, data = dat)))

dat %>% group_by(Sex) %>% nest() %>% 
  dplyr::mutate(
    t = map(data, function(d) {
      x = tidy(t.test(d$value ~ d$Condition))
      # x = tidy(kruskal.test(value ~ Condition, d))
    }
    )
  ) %>% dplyr::select(-data) %>% unnest()

dd = data.frame(Sex = c('Female', 'Male'), Sig = c('***', '*'), y = c(20, 24))
f2_1 = ggplot(dat, aes(x = Condition, y = ORG_BodyweightChange, fill = SexCondition)) +
  geom_hline(yintercept = 0, color = 'grey', lty = 2) +
  geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
  geom_point(pch = 21, size = 2, 
             position = position_jitterdodge(jitter.width = 0, 
                                             jitter.height = 0, dodge.width = 0.5)) + 
  scale_fill_manual(values = obj$four_color) + 
  scale_color_manual(values = obj$four_color) + 
  scale_y_continuous(breaks = seq(-30, 30, 10)) +
  labs(x = '', y = 'Bodyweight change') +
  facet_grid(.~Sex) + theme(legend.position = 'none') + 
  geom_text(data = dd, inherit.aes = F, aes(x = 1.5, y = y, label = Sig), fontface = 'bold', color = 'grey30')
f2_1

d2 = dat %>% group_by(Sex, Condition, SexCondition) %>% 
  dplyr::summarize(
    Md = median(ORG_FurScore_Cumulative, na.rm = T),
    MAD = stats::mad(ORG_FurScore_Cumulative, na.rm = T)
  )

dd = data.frame(Sex = c('Female', 'Male'), Sig = c('***', '***'), y = c(4, 6))

f2_2 = ggplot(dat, aes(x = Condition, y = ORG_FurScore_Cumulative, fill = SexCondition)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 1, position = position_dodge(0.7)) +
  scale_fill_manual(values = obj$four_color) + 
  scale_color_manual(values = obj$four_color) + 
  labs(x = '', y = 'Coat state') +
  facet_grid(.~Sex) + theme(legend.position = 'none') +
  geom_text(data = dd, inherit.aes = F, aes(x = 1.5, y = y, label = Sig), fontface = 'bold', color = 'grey30')
f2_2

dat$value = dat[['ORG_FurScore_Cumulative']]
dat %>% group_by(Sex, Condition) %>% 
  dplyr::summarize(
    n = n(), 
    n_groups = length(unique(GroupNumber)),
    M = mean(value, na.rm = T), 
    sd = sd(value, na.rm = T)
  )

dat %>% group_by(Sex, Condition) %>% nest() %>% 
  mutate(SW = map(data, function(d) shapiro.test(d$value)), 
         tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()

knitr::kable(car::leveneTest(value ~ Sex * Condition, dat))

dat = dat %>% 
  ungroup() %>% 
  # dplyr::mutate(value_trans = sqrt(value + 1)) %>%
  dplyr::mutate(value_trans = rank(value)) 

ggplot(dat, aes(x = Condition, y = value_trans, fill = SexCondition)) +
  geom_dotplot(binaxis='y', stackdir='center', dotsize = 1, position = position_dodge(0.7)) +
  scale_fill_manual(values = obj$four_color) + 
  scale_color_manual(values = obj$four_color) + 
  labs(x = '', y = 'Coat state') +
  facet_grid(.~Sex) + theme(legend.position = 'none')

dat %>% group_by(Sex, Condition) %>% nest() %>% 
  mutate(SW = map(data, function(d) shapiro.test(d$value_trans)), 
         tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()

tidy(kruskal.test(value ~ Condition, dat))

tidy(kruskal.test(value ~ Sex, dat))

dat %>% group_by(Sex) %>% 
  nest() %>% 
  dplyr::mutate(
    t = map(data, function(d) {
      tidy(kruskal.test(value ~ Condition, d))
    })
  ) %>% dplyr::select(-data) %>% unnest()

dd = data.frame(Sex = c('Female', 'Male'), Sig = c('', '**'), y = c(0, 0))

f2_3 = ggplot(dat, aes(x = Condition, y = ORG_Adrenals_adj, fill = SexCondition)) +
  geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
  geom_point(pch = 21, size = 2, 
             position = position_jitterdodge(jitter.width = 0, 
                                             jitter.height = 0, dodge.width = 0.5)) + 
  scale_fill_manual(values = obj$four_color) + 
  scale_color_manual(values = obj$four_color) + 
  labs(x = '', y = 'Adrenal weight') +
  facet_grid(.~Sex) + theme(legend.position = 'none') +
  geom_text(data = dd, inherit.aes = F, aes(x = 1.5, y = y, label = Sig), fontface = 'bold', color = 'grey30')
f2_3

dat$value = dat[['ORG_Adrenals_adj']]
dat %>% group_by(Sex, Condition) %>% 
  dplyr::summarize(
    n = n(), 
    n_groups = length(unique(GroupNumber)),
    M = mean(value, na.rm = T), 
    sd = sd(value, na.rm = T)
  )

dat %>% group_by(Sex, Condition) %>% nest() %>% 
  mutate(SW = map(data, function(d) shapiro.test(d$value)), 
         tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()

knitr::kable(car::leveneTest(value ~ Sex * Condition, dat))

as.data.frame(tidy(aov(value ~ Condition * Sex, data = dat)))

as.data.frame(tidy(aov(value ~ Condition + Sex, data = dat)))

dat %>% group_by(Sex) %>% nest() %>% 
  dplyr::mutate(
    t = map(data, function(d) {
      x = tidy(t.test(d$value ~ d$Condition))
      # x = tidy(kruskal.test(value ~ Condition, d))
    }
    )
  ) %>% dplyr::select(-data) %>% unnest()

f2_4 = ggplot(dat, aes(x = Condition, y = ORG_Cort, fill = SexCondition)) +
  # geom_hline(yintercept = 0, color = 'grey', lty = 2) +
  geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
  geom_point(pch = 21, size = 2, 
             position = position_jitterdodge(jitter.width = 0, 
                                             jitter.height = 0, dodge.width = 0.5)) + 
  scale_fill_manual(values = obj$four_color) + 
  scale_color_manual(values = obj$four_color) + 
  labs(x = '', y = 'Corticosterone') +
  facet_grid(.~Sex) + theme(legend.position = 'none')
f2_4

dat$value = dat[['ORG_Cort']]
dat %>% group_by(Sex, Condition) %>% 
  dplyr::summarize(
    n = n(), 
    n_groups = length(unique(GroupNumber)),
    M = mean(value, na.rm = T), 
    sd = sd(value, na.rm = T)
  )

dat %>% group_by(Sex, Condition) %>% nest() %>% 
  mutate(SW = map(data, function(d) shapiro.test(d$value)), 
         tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()

knitr::kable(car::leveneTest(value ~ Sex * Condition, dat))

knitr::kable(car::leveneTest(value ~ Condition, dat[dat$Sex == 'Male',]))

knitr::kable(car::leveneTest(value ~ Condition, dat[dat$Sex == 'Female',]))

tidy(kruskal.test(value ~ Condition, dat))

tidy(kruskal.test(value ~ Sex, dat))

dat %>% group_by(Sex) %>% 
  nest() %>% 
  dplyr::mutate(
    t = map(data, function(d) {
      tidy(kruskal.test(value ~ Condition, d))
    })
  ) %>% dplyr::select(-data) %>% unnest()

Correlation plot

dat = obj$tests_combined_adj
props = colnames(dat)[regexpr('_', colnames(dat)) > 0]

props_pretty = c(
  'EPM | Center time', 
  'EPM | Closed arm distance', 
  'EPM | Closed arm entries', 
  'EPM | Closed arm mean visit',  
  'EPM | Closed arm time', 
  'EPM | Distance', 
  'EPM | Open arm distance from center', 
  'EPM | Open arm distance', 
  'EPM | Open arm entries', 
  'EPM | Open arm mean visit', 
  'EPM | Open arm time', 
  'EPM | Closed arm preference', 
  
  'NBT | Intact material', 
  'NBT | Nest score', 
  
  'OFT | Mean distance from center', 
  'OFT | Distance',
  'OFT | Inner zone mean speed', 
  'OFT | Inner zone distance', 
  'OFT | Inner zone visits',
  'OFT | Inner zone mean visit duration', 
  
  'OFT | Max. speed', 
  'OFT | Meander',
  'OFT | Outer zone mean speed', 
  'OFT | Outer zone distance', 
  'OFT | Outer zone mean visit duration',
  
  'OFT | Outer zone preference', 
  'OFT | Immobility time',
  
  'PHYS | Adrenal weight',
  'PHYS | Bodyweight change', 
  'PHYS | Corticosterone', 
  'PHYS | Coat state', 
  
  'SPT | Median preference', 
  'SPT | Total preference', 
  
  'SPL | Grooming time',
  'SPL | Latency to groom', 
  'TST | Immobility episodes',
  'TST | Immobility time'
)

as.data.frame(cbind(props, props_pretty))

pcaDat = dat[,colnames(dat) %in% props]
pheno = dat[,!colnames(dat) %in% props]

h = as.matrix(pcaDat)
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)

hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')

f2_5 = ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = 'white') +
  scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
  scale_y_discrete(position = 'right') +
  labs(x = '', y = '') + 
  coord_fixed() +
  theme(
    axis.text.x = element_blank(), panel.border = element_blank(),
    plot.margin = margin(15,0,0,15, "pt"), 
    legend.position = 'bottom'
  )
f2_5

# MEs = pcaDat
# plot = pheatmap::pheatmap(cor(MEs), border_color = 'white', kmeans_k = 3,
# fontsize = 7,  treeheight_col = 0) 
# plot

Dimensions

dim(pcaDat)

## [1] 73 37

Samples

pheno %>% group_by(Sex, Condition) %>% dplyr::summarise(N = n())

Figure panel

row1 = cowplot::plot_grid(plotlist = list(f2_1, f2_2, f2_3, f2_4), ncol = 2, align = 'hv',
                          labels = c('b', 'c', 'd', 'e'), label_fontface = 'plain')
row2 = cowplot::plot_grid(plotlist = list(f2_5), 
                          labels = c('f'), label_fontface = 'plain')
cowplot::plot_grid(plotlist = list(row1, row2), nrow = 2, 
                   rel_heights = c(1, 1.2), axis = 'tb') 

if(exportFigures) {
  ggsave('../manuscript/figures/2_Panel.pdf', width = 7, height = 10, dpi = 1200)
  ggsave('../manuscript/figures/2_Panel.png', width = 7, height = 10, dpi = 1200)
}

PCA

• Adjusted for batch effects

dat = obj$tests_combined_adj

pcaDat = dat[,colnames(dat) %in% props]
pheno = dat[,!colnames(dat) %in% props]

# pca = prcomp(pcaDat, center = T, scale. = F)
pca = prcomp(pcaDat, center = T, scale. = T)
PCs = colnames(as.data.frame(pca$x[,1:6]))
p = cbind(as.data.frame(pheno), as.data.frame(pca$x[,1:6]))

p$Sex = factor(p$Sex, levels = c('Female', 'Male'))
p$SexCondition = paste0(p$Sex, '_', dat$Condition)
p$SexCondition = factor(p$SexCondition, 
                        levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))
obj$analyses$PCA_OutcomeVars_BatchAdj = p

Variance explained

eigs <- pca$sdev^2
e = data.frame(
  PC = colnames(pca$x),
  Proportion = eigs/sum(eigs),
  Cumulative = cumsum(eigs)/sum(eigs))
e = e[1:5,] %>% 
  mutate(PC = factor(PC, levels = colnames(pca$x)))
f3_1 = ggplot(e, aes(x = PC, y = Proportion*100)) + 
  geom_line(aes(group = 1), color = 'grey') +
  geom_point(size = 2, pch = 21, fill = obj$colors[2]) +
  scale_y_continuous(limits = c(5, 25), breaks = seq(0, 25, 5)) +
  labs(x = '', y = 'Variance explained (%)')
f3_1

as_tibble(e)

PCs vs sex & condition

Check Normality

idVars = c(colnames(pheno), 'SexCondition')
PCs = c('PC1', 'PC2', 'PC3', 'PC4')

pp <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  dplyr::select(idVars, PCs) %>%
  gather(PC, value, -idVars) %>%
  group_by(PC, Sex, Condition) %>% nest() %>%
  mutate(SW = map(data, function(d) shapiro.test(d$value)), 
         tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()

Normality is not violated in any groups

Check heteroskedasticity

pp <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  dplyr::select(idVars, PCs) %>%
  gather(PC, value, -idVars) %>%
  group_by(PC) %>% nest() %>%
  dplyr::mutate(
    t = map(data, function(d) {
      car::leveneTest(value ~ Sex * Condition, d)
    })
  ) %>% dplyr::select(-data) %>% unnest()
pp

All variances are homogenous

idVars = c(colnames(pheno), 'SexCondition')
PCs = c('PC1', 'PC2', 'PC3', 'PC4')

pp <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  dplyr::select(idVars, PCs) %>%
  gather(PC, value, -idVars) %>%
  group_by(PC) %>% nest() %>%
  mutate(
    m = map2(data, PC, function(d, p) {
      summary(aov(value ~ Condition * Sex, d))
    })
  ) %>% dplyr::select(-data) 

plots = obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  dplyr::select(idVars, PCs) %>%
  gather(PC, value, -idVars) %>%
  group_by(PC) %>% nest() %>% mutate(
    plot = map2(data, PC, function(d, p) {
      ggplot(d, aes(x = Condition, y = value, fill = SexCondition)) +
        geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition),
                     position = position_dodge(0.6)) +
        geom_point(pch = 21, position = position_dodge(0.6), size = 2.5) +
        labs(x = '', y = p) + facet_grid(~Sex) +
        scale_fill_manual(values = obj$four_color) +
        scale_color_manual(values = obj$four_color) +
        theme(legend.position = 'none')
    })
  )

PC1:

pp$m[[1]]

##               Df Sum Sq Mean Sq F value Pr(>F)
## Condition      1    3.3   3.304   0.411  0.523
## Sex            1    0.1   0.140   0.017  0.895
## Condition:Sex  1    8.4   8.401   1.046  0.310
## Residuals     69  554.3   8.033

PC2:

pp$m[[2]]

##               Df Sum Sq Mean Sq F value         Pr(>F)    
## Condition      1   3.34    3.34   1.233          0.271    
## Sex            1 139.50  139.50  51.449 0.000000000646 ***
## Condition:Sex  1   6.04    6.04   2.226          0.140    
## Residuals     69 187.09    2.71                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PC3:

pp$m[[3]]

##               Df Sum Sq Mean Sq F value  Pr(>F)   
## Condition      1  28.70  28.699   7.897 0.00644 **
## Sex            1   1.10   1.097   0.302 0.58449   
## Condition:Sex  1   0.87   0.874   0.241 0.62530   
## Residuals     69 250.74   3.634                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PC4:

pp$m[[4]]

##               Df Sum Sq Mean Sq F value  Pr(>F)   
## Condition      1   0.42   0.417   0.159 0.69176   
## Sex            1  27.24  27.237  10.344 0.00198 **
## Condition:Sex  1  30.08  30.076  11.423 0.00120 **
## Residuals     69 181.68   2.633                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

grid.arrange(grobs = plots$plot, ncol = 2)

f3_2 = plots$plot[plots$PC == 'PC1'][[1]] + ylab('PC1 score')

PC1

baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                             idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(as.character(MouseID)))
pp = left_join(baseline, obj$analyses$PCA_OutcomeVars_BatchAdj)
pp$Condition = factor(pp$Condition, levels = c('Control', 'CMS'))
f3_3 = ggplot(pp, aes(x = DS, y = PC1)) +
  geom_smooth(method = 'lm', se = F, lty = 2, aes(col = SexCondition)) +
  geom_point(pch = 21, size = 2, aes(fill = SexCondition)) +
  labs(y = "PC1 score", x = "Baseline David's Score") +
  facet_grid(Condition~Sex) +
  scale_fill_manual(values = obj$four_color) +
  scale_color_manual(values = obj$four_color) +
  theme(legend.position = 'none')
f3_3

Stats

summary(lm(PC1 ~ DS * Sex, pp[pp$Condition == 'CMS',]))

## 
## Call:
## lm(formula = PC1 ~ DS * Sex, data = pp[pp$Condition == "CMS", 
##     ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.9194 -1.7251 -0.1317  1.6727  5.5421 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   -4.421      2.045  -2.162  0.03754 * 
## DS             2.856      1.287   2.220  0.03301 * 
## SexMale        9.246      2.591   3.568  0.00107 **
## DS:SexMale    -5.834      1.637  -3.563  0.00108 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.554 on 35 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.2857, Adjusted R-squared:  0.2245 
## F-statistic: 4.667 on 3 and 35 DF,  p-value: 0.007594

summary(aov(PC1 ~ DS * Sex, pp[pp$Condition == 'CMS',]))

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## DS           1   6.22    6.22   0.954 0.33533   
## Sex          1   2.30    2.30   0.353 0.55639   
## DS:Sex       1  82.78   82.78  12.695 0.00108 **
## Residuals   35 228.22    6.52                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 9 observations deleted due to missingness

p = obj$analyses$PCA_OutcomeVars_BatchAdj

PC = 'PC1'
p$PC = p[[PC]]
cors = sapply(1:ncol(pcaDat), function(prop) { cor(p$PC, pcaDat[,prop], method = 'spearman') })
pvals = sapply(1:ncol(pcaDat), function(prop) { cor.test(p$PC, pcaDat[,prop], method = 'spearman')$p.value })

cors = data.frame(feature = colnames(pcaDat), rho = cors, pval = pvals)
cors$feature =  factor(cors$feature, levels = cors$feature[order(abs(cors$rho), decreasing = F)])
cors = cors[order(abs(cors$rho), decreasing = T), ]
cors$Sign = ifelse(cors$rho < 0, yes = 'neg', no = 'pos')
cors$Sign = factor(cors$Sign, levels = c('pos', 'neg'))

# cors$Sig = cors$pval <= 0.05
# cors$Sig = p.adjust(cors$pval, method = 'BH') <= 0.1

# ? Bonferroni under 0.1 - is that good enough ? 

cors$Sig = p.adjust(cors$pval, method = 'bonferroni') <= 0.1
sum(cors$Sig)

## [1] 16

cors$feature_pretty = plyr::mapvalues(cors$feature, props, props_pretty)

hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')

f3_4 = ggplot(cors, aes(y = abs(rho), x = feature_pretty, col = Sign, fill = Sign)) +
  geom_bar(stat = 'identity', width = 0.01) +  
  geom_point(size = 2, pch = 21) + 
  geom_point(data = subset(cors, Sig), size = 2.5, pch = 21, color = 'black') +
  scale_fill_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = '') +
  scale_color_manual(values = c(hm.palette(9)[8], hm.palette(9)[1]), name = '') +
  # xlab('') + ylab() + coord_flip() + obj$theme +
  xlab('') + ylab("Correlation") + coord_flip() + obj$theme +
  theme(plot.title = element_text(hjust = 0.5), legend.position = 'top')  
f3_4

dat = obj$tests_combined_adj
idVars = c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber', 
           'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props = colnames(dat)[!colnames(dat) %in% idVars]
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                             idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(as.character(MouseID))) %>% 
  rename(Dominance = DS)
d = left_join(dat, baseline)

idVars = c(idVars, 'Dominance')
plots = d  %>% 
  dplyr::select(idVars, props) %>% 
  gather(prop, value, -idVars) %>% 
  group_by(prop) 
plots$prop = plyr::mapvalues(plots$prop, props, props_pretty)
plots$Sex = factor(plots$Sex, levels = c('Female', 'Male'))
plots$Condition = factor(plots$Condition, levels = c('Control', 'CMS'))
plots$SexCondition = paste0(plots$Sex, '_', dat$Condition)
plots$SexCondition = factor(plots$SexCondition, 
                            levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))

selected_props = c('OFT | Distance', 'EPM | Closed arm time')
# selected_props = c('OFT | Inner zone distance', 'EPM | Open arm distance')

res = lapply(selected_props, function(i) {
  pp = subset(plots, prop == i)
  l = c(min(pp$value), max(pp$value))
  p1 = ggplot(pp, aes(x = Condition, y = value, fill = SexCondition)) +
    geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
    geom_point(pch = 21, size = 2, 
               position = position_jitterdodge(jitter.width = 0, 
                                               jitter.height = 0, dodge.width = 0.5)) + 
    scale_fill_manual(values = obj$four_color) + 
    scale_color_manual(values = obj$four_color) + 
    facet_grid(Sex~.) + labs(x = '', y = i) +
    scale_y_continuous(limits = l) +
    theme(legend.position = 'none') + 
    theme(strip.text.y = element_blank())
  
  p2 = ggplot(pp[pp$Condition == 'CMS',], aes(x = Dominance, y = value, fill = Sex, color = Sex)) +
    geom_smooth(method = 'lm', se = F, lty = 2) +
    geom_point(pch = 21, size = 2, color = 'grey20') +
    scale_fill_manual(values = obj$four_color[3:4]) +
    scale_color_manual(values = obj$four_color[3:4]) +
    scale_y_continuous(limits = l) +
    labs(y = '', x = "Baseline David's Score") +
    facet_grid(Sex~Condition) +
    theme(legend.position = 'none', 
          axis.text.y = element_blank(),
          plot.margin = margin(0,0,0,-10, "pt")
    )
  
  cowplot::plot_grid(plotlist = list(p1, p2), axis = 'tb',
                     align = 'h', rel_widths = c(1,1.5))
})
f3_5 = cowplot::plot_grid(plotlist = res, ncol = 2, align = 'v', rel_widths = c(1,1),
                          labels = c('e', 'f'), label_fontface = 'plain')
f3_5

Figure panel

row1 = cowplot::plot_grid(plotlist = list(f3_1, f3_2), ncol = 2, axis = 't', align = 'h',
                          labels = c('a', 'b'), label_fontface = 'plain', rel_widths = c(0.7,1))

row2 = cowplot::plot_grid(plotlist = list(f3_3), ncol = 1, align = 'hv',
                          labels = c('c'), label_fontface = 'plain')

col1 = cowplot::plot_grid(plotlist = list(row1, row2), ncol = 1, axis = 't', align = 'h',
                          rel_widths = c(0.7,1), rel_heights = c(1, 2.1))

col2 = cowplot::plot_grid(plotlist = list(f3_4), ncol = 1, align = 'hv',
                          labels = c('d'), label_fontface = 'plain')

l1 = cowplot::plot_grid(plotlist = list(col1, col2), align = 'hv', ncol = 2, rel_widths = c(1, 0.9),
                        label_fontface = 'plain')

cowplot::plot_grid(plotlist = list(l1, f3_5), align = 'hv', ncol = 1, rel_heights = c(1, 0.6))

if(exportFigures) {
  ggsave('../manuscript/figures/3_Panel.pdf', width = 9, height = 11, dpi = 1200)
  ggsave('../manuscript/figures/3_Panel.png', width = 9, height = 11, dpi = 1200)
}

PC1 versus ranks

idVars = c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber', 
           'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props = colnames(dat)[!colnames(dat) %in% idVars]
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                             idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(as.character(MouseID))) %>% 
  rename(Dominance = DS) %>% group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-Dominance), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )

baseline = left_join(baseline, obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])
p <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  dplyr::select(RealMouseNumber, PC1)

d = left_join(baseline, p)
d$Sex = factor(d$Sex, levels = c('Female', 'Male'))
d$Condition = factor(d$Condition, levels = c('Control', 'CMS'))
d$SexCondition = paste0(d$Sex, '_', dat$Condition)
d$SexCondition = factor(d$SexCondition, 
                        levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))
d$Rank_3 = plyr::mapvalues(d$Rank, c('alpha', 'beta', 'gamma', 'delta'), 
                           c('alpha', 'beta / gamma', 'beta / gamma', 'delta'))
d$Rank_3 = factor(d$Rank_3, levels = c('alpha', 'beta / gamma', 'delta'))
d$Rank_3_ordinal = as.numeric(as.character(plyr::mapvalues(d$Rank_3,  
                                                           c('alpha', 'beta / gamma', 'delta'), 
                                                           c(1, 0, -1))))
pp = d
pp$value = pp$PC1
l = c(min(pp$value), max(pp$value))

xx = data.frame(
  Sex = c('Female', 'Female', 'Male', 'Male'), 
  x = c(1.5, 2.5, 1.5, 2.5),
  y = c(5, 0, 4.5, 5),
  Sig =  c('*', '', '#', '*'), 
  Line =  c('_', '', '_', '_')
  )
f4_1 = ggplot(pp[pp$Condition == 'CMS',], aes(x = Rank_3, y = value, fill = Rank_3)) +
  geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = Rank_3)) +
  geom_point(pch = 21, size = 2, 
             position = position_jitterdodge(jitter.width = 0, 
                                             jitter.height = 0, dodge.width = 0.5)) + 
  geom_text(data = xx[-3,], inherit.aes = F, aes(x = x, y = y, label = Sig), fontface = 'bold', size = 4,
            color = 'grey30') +
  geom_text(data = xx[3,], inherit.aes = F, aes(x = x, y = y, label = Sig), size = 3, fontface = 'bold', color = 'grey30') +
  
  geom_text(data = xx, inherit.aes = F, aes(x = x, y = y, label = Line), fontface = 'bold', color = 'grey30') +
  
  scale_color_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4])) +
  scale_fill_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4]), 
                    labels = c('\u03B1', paste('\u03B2', '\u03B3', sep = '/'), '\u03B4')) +
  scale_x_discrete(labels = c('\u03B1', paste('\u03B2', '\u03B3', sep = '/'), '\u03B4')) +
  scale_y_continuous(limits = l) +
  facet_grid(Sex~Condition) + labs(x = '', y = 'PC1 score') +
  theme(legend.position = 'none',
        axis.title.x = element_text(family="Helvetica")
  ) 
f4_1

dat = subset(pp, Condition == 'CMS')
dat$value = dat[['PC1']]
dat %>% group_by(Sex, Rank_3) %>% 
  filter(!is.na(value)) %>%
  dplyr::summarize(
    n = n(), 
    n_groups = length(unique(GroupNumber)),
    M = mean(value, na.rm = T), 
    sd = sd(value, na.rm = T)
  )

dat %>% group_by(Sex, Rank_3) %>% nest() %>% 
  mutate(SW = map(data, function(d) shapiro.test(d$value)), 
         tSW = map(SW, broom::tidy)) %>% select(-data, -SW) %>% unnest()

knitr::kable(car::leveneTest(value ~ Sex * Rank_3, dat))

dds = dat %>% 
  filter(Sex == 'Male')
as.data.frame(summary(aov(value ~ Rank_3_ordinal, dds))[[1]])

dds = dat %>% 
  filter(Sex == 'Male' & Rank_3_ordinal %in% c(0, 1))
dds$Rank_3 = droplevels(dds$Rank_3)
as.data.frame(summary(aov(value ~ Rank_3, dds))[[1]])

dds = dat %>% 
  filter(Sex == 'Male' & Rank_3_ordinal %in% c(0, -1))
dds$Rank_3 = droplevels(dds$Rank_3)
as.data.frame(summary(aov(value ~ Rank_3, dds))[[1]])

dds = dat %>% 
  filter(Sex == 'Female')
as.data.frame(summary(aov(value ~ Rank_3_ordinal, dds))[[1]])

dds = dat %>% 
  filter(Sex == 'Female' & Rank_3_ordinal %in% c(0, 1))
dds$Rank_3 = droplevels(dds$Rank_3)
as.data.frame(summary(aov(value ~ Rank_3, dds))[[1]])

dds = dat %>% 
  filter(Sex == 'Female' & Rank_3_ordinal %in% c(0, -1))
dds$Rank_3 = droplevels(dds$Rank_3)
as.data.frame(summary(aov(value ~ Rank_3, dds))[[1]])

Current

Select top PC1 behaviors

• All features with with a significant adjusted PC1 correlation

# selected_props = as.character(cors$feature[abs(cors$rho) >= 0.5])
selected_props = as.character(cors$feature[cors$Sig == T])

dat = obj$tests_combined_adj
props = colnames(dat)[regexpr('_', colnames(dat)) > 0]

pcaDat = dat[,colnames(dat) %in% props]
pheno = dat[,!colnames(dat) %in% props]

pcaDat = pcaDat[,colnames(pcaDat) %in% selected_props]

h = as.matrix(pcaDat)
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)

hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')

ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = 'white') +
  scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
  scale_y_discrete(position = 'right') +
  labs(x = '', y = '') + 
  coord_fixed() +
  theme(
    axis.text.x = element_blank(), panel.border = element_blank(),
    plot.margin = margin(15,0,0,15, "pt"), 
    legend.position = 'bottom'
  )

All Controls

h = as.matrix(pcaDat[pheno$Condition == 'Control',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = 'white') +
  scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
  scale_y_discrete(position = 'right') +
  labs(x = '', y = '') + 
  coord_fixed() +
  theme(
    axis.text.x = element_blank(), panel.border = element_blank(),
    plot.margin = margin(15,0,0,15, "pt"), 
    legend.position = 'bottom'
  )

All CMS

h = as.matrix(pcaDat[pheno$Condition == 'CMS',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = 'white') +
  scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
  scale_y_discrete(position = 'right') +
  labs(x = '', y = '') + 
  coord_fixed() +
  theme(
    axis.text.x = element_blank(), panel.border = element_blank(),
    plot.margin = margin(15,0,0,15, "pt"), 
    legend.position = 'bottom'
  )

h = as.matrix(pcaDat[pheno$Condition == 'CMS',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)

# plot(hc, hang = -1, cex = 0.6)
hc$labels = plyr::mapvalues(hc$labels, props, props_pretty)
hcd <- as.dendrogram(hc)
# Default plot
# plot(hcd, type = "rectangle", ylab = "Height")
ggdendro::ggdendrogram(hc, rotate = T)

# 
# dend <- as.dendrogram(hc)
# # Extract the data (for rectangular lines)
# # Type can be "rectangle" or "triangle"
# dend_data <- ggdendro::dendro_data(dend, type = "rectangle")
# ggplot(dend_data$segments) + 
#   geom_segment(aes(x = x, y = y, xend = xend, yend = yend))

Male Controls

h = as.matrix(pcaDat[pheno$Sex == 'Male' & pheno$Condition == 'Control',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = 'white') +
  scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
  scale_y_discrete(position = 'right') +
  labs(x = '', y = '') + 
  coord_fixed() +
  theme(
    axis.text.x = element_blank(), panel.border = element_blank(),
    plot.margin = margin(15,0,0,15, "pt"), 
    legend.position = 'bottom'
  )

Male CMS

h = as.matrix(pcaDat[pheno$Sex == 'Male' & pheno$Condition == 'CMS',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = 'white') +
  scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
  scale_y_discrete(position = 'right') +
  labs(x = '', y = '') + 
  coord_fixed() +
  theme(
    axis.text.x = element_blank(), panel.border = element_blank(),
    plot.margin = margin(15,0,0,15, "pt"), 
    legend.position = 'bottom'
  )

Female Control

h = as.matrix(pcaDat[pheno$Sex == 'Female' & pheno$Condition == 'Control',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = 'white') +
  scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
  scale_y_discrete(position = 'right') +
  labs(x = '', y = '') + 
  coord_fixed() +
  theme(
    axis.text.x = element_blank(), panel.border = element_blank(),
    plot.margin = margin(15,0,0,15, "pt"), 
    legend.position = 'bottom'
  )

Female CMS

h = as.matrix(pcaDat[pheno$Sex == 'Female' & pheno$Condition == 'CMS',])
cormat = round(cor(h),2)
hc = hclust(as.dist(1-cormat)/2)
cormat.ord = cormat[hc$order, hc$order]
melted_cormat = reshape2::melt(cormat.ord, na.rm = T)
melted_cormat$Var1 = plyr::mapvalues(melted_cormat$Var1, props, props_pretty)
melted_cormat$Var2 = plyr::mapvalues(melted_cormat$Var2, props, props_pretty)
hm.palette = colorRampPalette(rev(RColorBrewer::brewer.pal(9, 'Spectral')), space='Lab')
ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = 'white') +
  scale_fill_gradientn(colours = hm.palette(100), name = 'Correlation') +
  scale_y_discrete(position = 'right') +
  labs(x = '', y = '') + 
  coord_fixed() +
  theme(
    axis.text.x = element_blank(), panel.border = element_blank(),
    plot.margin = margin(15,0,0,15, "pt"), 
    legend.position = 'bottom'
  )